//============================================================================
// Name        : Electrostatics.cpp
// Author      : 
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//============================================================================

#include <iostream>
#include <fstream>
#include <sstream>
#include "Rho.h"
#include "RhoDip.h"
#include "typedefs.h"
#include <string>
using namespace std;
enum {Slt,Sol,Tot};
#include "Phi.h"
#include "Dielectric.h"
#include "DiffCoeffs_Spec.hpp"
#include "Grid_SpecTemplates.hpp"
#include <vector>
#include "EpsofK.h"
#include "simpson.hpp"

int main(int argc, char **argv) {
	string filename("Rhos.bin");

	int nFilter=1;;
	bool bFilter=FALSE;
	int n=1;
	cout << " Out of " << endl;
	std::stringstream readit;

	for(;;){
		if(n>=argc) break;
		string opt;
		opt.assign(argv[n]);

		if(opt[0] == '-') {
			switch(opt[1]){
			case 'h':
				cout<< "Usage: " << argv[0] << " -i infile -f nFilter " << endl;
				exit(1);
			case 'i':
				filename.assign(argv[++n]);
				break;
			case 'f':
				bFilter=TRUE;
				readit<< argv[++n];
				readit >> nFilter;
				break;
			default:
				cout << "option not found abort " << endl;
				exit(1);
			}
		}
		n++;
	}

    Rho Rhos[DIM],Rho_Tot;

	RhoDip Pol;
	ifstream fin;
	fin.open(filename.c_str(),std::ios_base::in|std::ios_base::binary);

	fin >> Rhos[Slt]>> Rhos[Sol] >> Rhos[Tot] >> Pol;

	matrix & co=Rhos[Sol].getCO();
	rvec x={0.5*co[0][0],0.5*co[1][1],0.5*co[2][2]};
	FILE * fg4=fopen("MyRhos.xvg","w+");

	double x2[DIM];
	for(int i=0;i<DIM;i++) x2[i]=x[i];
	Rhos[Sol].Rdf(fg4,x2,"Rhos Sol");
	exit(1);

//	Rhos[Sol].Symmetrify(x);
	EpsofK Eps0;
	Eps0(Rhos[Slt],Rhos[Sol]);
	FILE * fg3=fopen("MyElectrostatics.xvg","w+");
	Eps0.RdfK(fg3,"Eps ");

	exit(1);



	if(bFilter) {
		Grid<1>::MakeFilter(nFilter,nFilter,nFilter);
		Grid<DIM>::MakeFilter(nFilter,nFilter,nFilter);
	}
    double Myx[DIM]={Grid<1>::getCO()[0][0]*0.5,Grid<1>::getCO()[1][1]*0.5,Grid<1>::getCO()[2][2]*0.5 };
	FILE * fg1=fopen("MyFieldsn.xvg","w+");
	Phi MyPhi[3];
	Phi Phi_tot;
    Rho Rho_p;
    Phi Phi_p;
    Rho_Tot.Allocate();
	Rhos[Sol].Filter();
	Rhos[Slt].Filter();
	Pol.Filter();
	Rhos[Tot]=Rhos[Sol]+Rhos[Slt];
    Rho_p=Pol.Div();
    Rho_p=-Rho_p;
    Rho_Tot=Rho_p+Rhos[Slt];

	int nfx,nfy,nfz;
	int nx0=static_cast<int>(Grid<1>::getnnx());
	int ny0=static_cast<int>(Grid<1>::getnny());
	int nz0=static_cast<int>(Grid<1>::getnnz());
    Grid<DIM> E_Slt, E_p,E_Tot,E_Sol, E;

	nfx=(nx0 % 2 == 0)? nx0/2: nx0/2+1;
	nfy=(ny0 % 2 == 0)? ny0/2: ny0/2+1;
	nfz=(nz0 % 2 == 0)? nz0/2: nz0/2+1;


//	matrix & co=Rhos[Sol].getCO();
//	rvec x={0.5*co[0][0],0.5*co[1][1],0.5*co[2][2]};
	Rhos[Sol].Symmetrify(x);

    MyPhi[Slt](Rhos[Slt]);
    MyPhi[Sol](Rhos[Sol]);
    MyPhi[Tot]=MyPhi[Slt]+MyPhi[Sol];
    std::vector<double> Phi_t = MyPhi[Tot].Rdf(fg1,Myx,"Phi_t");
    std::vector<double> Phi_v0 = MyPhi[Slt].Rdf(fg1,Myx,"Phi_0");
    cout << " & " << endl;
    double dx=Phi_t[0];
    for(int i=1;i<Phi_t.size();i++){
    	double x=(i-1)*dx;
    	double fx=(Phi_v0[i])?Phi_t[i]/Phi_v0[i]:0.0;
    	cout << x << " " << fx << endl;
    }
    cout << " & " << endl;
    exit(1);
    Phi_p(Rho_p);
    Phi_tot=Phi_p+MyPhi[Slt];

    E_Tot=MyPhi[Tot].Field();
    E_Slt=MyPhi[Slt].Field();
    E_Sol=MyPhi[Sol].Field();
    E_p=Phi_p.Field();

    E=E_Slt+E_p;

    std::vector<double> Rho_Slt = Rhos[Slt].Rdf(fg1,Myx,"Rho_Slt");
    std::vector<double> Rho_Sol = Rhos[Sol].Rdf(fg1,Myx,"Rho_Sol");
    std::vector<double> Rho_p0 = Rho_p.Rdf(fg1,Myx,"Rho_p0");
//    std::vector<double> Rho_To = Rhos[Tot].Rdf(fg1,Myx,"Rho_Tot");
    std::vector<double> Pol_0 = Pol.Rdf(fg1,Myx,"Pol");

//    std::vector<double> Phi_vtot = MyPhi[Tot].Rdf(fg1,Myx,"Phi_t");
//    std::vector<double> Phi_vsol = MyPhi[Sol].Rdf(fg1,Myx,"Phi_s");
//    std::vector<double> Phi_v0 = MyPhi[Slt].Rdf(fg1,Myx,"Phi_0");

//    std::vector<double> EE = E.Rdf(fg1,Myx,"E");
//    std::vector<double> EE_Tot = E_Tot.Rdf(fg1,Myx,"E_Tot");
//    std::vector<double> EE_Slt = E_Slt.Rdf(fg1,Myx,"E_Slt");
    Dielectric test(nx0,ny0,nz0);
	FILE * fg2=fopen("MyEps.xvg","w+");
//    test.GetEps(fg2,"E",E,E_Slt);
//    test.GetEps(fg2,"E_tot",E_Tot,E_Slt);
    test.GetEps(fg2,"Rhos_Tot Rhos_slt",Rhos[Tot],Rhos[Slt]);

/*


    for(unsigned int i=0;i<Grid<DIM>::getnnx();i++)
		cout << i << " " << Pol[ZZ][nfx][nfy][i] << " " << E_p[ZZ][nfx][nfy][i]<< " " <<
			E_Slt[ZZ][nfx][nfy][i] << " " << E_Sol[ZZ][nfx][nfy][i]<<
			" " << E_Tot[ZZ][nfx][nfy][i]<< " " << E[ZZ][nfx][nfy][i]<< endl;
*/

/*
    std::vector<double> V_p = Pol.Rdf(fg1,Myx,"Pol");
    std::vector<double> V_eSlt = E_Slt.Rdf(fg1,Myx,"E_0");
    std::vector<double> V_e = E.Rdf(fg1,Myx,"E");
    std::vector<double> V_eTot = E_Tot.Rdf(fg1,Myx,"E_Tot");
    std::vector<double> V_ePol = E_p.Rdf(fg1,Myx,"E_Pol");
    std::vector<double> V_eSol = E_Sol.Rdf(fg1,Myx,"E_Pol");
*/
/*
    double dx=V_eSlt[0];
    for(unsigned int i=1;i<V_eSlt.size();i++){
    	double eps1=1.0+(4.0*M_PI)*V_p[i]/V_eTot[i];
    	double eps2=1.0+(4.0*M_PI)*V_p[i]/V_e[i];
    	double eps3=V_eSlt[i]/V_eTot[i];
    	cout << (i-1)*dx << "  " << eps1 << " " << eps2 << " " << eps3 << endl;
    }

*/
    exit(1);
	return 0;
}
